#!/usr/bin/python
# GM_velocities.py
'''Create velocity distributions for the Milky Way.'''

import math
import numpy
from matplotlib import use as useBackend
useBackend('Agg')
import matplotlib.pyplot as ppt


def cheat111():
    '''Cheat to clear plot.'''
    ppt.subplot(311)
    ppt.plot([0,1],[1,2])
    ppt.subplot(131)
    ppt.plot([0,1],[1,2])
    ppt.subplot(111)


VariableNames = ['d', 'la', 'ba', 'phi', 'bulk_vels', 'sigma_vels']

def velocities(simpars, phypars):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    # Relevant distances
    # 1d
    d2 = d**2
    cosl = numpy.cos(la)
    sinl = numpy.sin(la)
    cosb = numpy.cos(ba)
    sinb = numpy.sin(ba)
    # 3d
    cosphi = numpy.cos(phi)
    sinphi = numpy.sin(phi)
    alpha = numpy.subtract(phi, la[numpy.newaxis,:,numpy.newaxis])
    cosalpha = numpy.cos(alpha)
    sinalpha = numpy.sin(alpha)
#    del phi
    del alpha
    # I keep some "vestigial" terms in comments that are presently useless but also harmless, which (technically) may be of future use.
    # Bulk motion first.
    # Rg, phi, z to x, y, z.
    (vabsRg, vabsphi, vabsz) = bulk_vels
# ---> HERE is where I should insert a more detailed description of rotation velocity when I get around to that.  Just replace vabsphi (currently a lone float) with a 3d array by the same name, drawing on R (which is premade and easily imported and all).  The quantities vabsphi gets multiplied by, sinphi and cosphi, are alread 3d so no extra sweat is required.
    # Sign here "flipped" from typical convention b/c galaxy rotates clockwise about north galactic pole, "left-handed"
    vabsx = sinphi*vabsphi #+ any R term # vestigial
    vabsy = -cosphi*vabsphi #+ any R term # vestigial
    #vabsz = vabsz # any z term # vestigial
    # Here I remove the sun's own absolute velocity, which is in the +y/-phi direction.
    (vx0, vy0, vz0) = bulk_vels
    vrelx = vabsx# - vx0 # vestigial
    vrely = vabsy - vy0
    vrelz = vabsz# - vz0 # vestigial
    # x, y, z direct to d, l, b; all coords give 3d result
    #vreld = cosb*cosl*vrelx + cosb*sinl*vrely + sinb*vrelz # vestigial
    vrell = -numpy.multiply(sinl[numpy.newaxis,:,numpy.newaxis], vrelx) + numpy.multiply(cosl[numpy.newaxis,:,numpy.newaxis], vrely)
    vrelappl = numpy.divide(vrell, d[:,numpy.newaxis,numpy.newaxis])
    #vrelb = -sinb*cosl*vrelx - sinb*sinl*vrely + cosb*vrelz # vestigial
    # Now dispersed motion.
    # Direct without intermediaries to final coord system!
    (sigmaRg, sigmaphi, sigmaz) = sigma_vels
    #sigmad2 = (cosb*cosalpha*sigmaRg)**2 + (-cosb*sinalpha*sigmaphi)**2 + (sinb*sigmaz)**2 # vestigial
    sigmal2 = (sinalpha*sigmaRg)**2 + (cosalpha*sigmaphi)**2
    sigmab2 = ( numpy.multiply(-sinb[numpy.newaxis,numpy.newaxis,:], cosalpha*sigmaRg) )**2 + ( numpy.multiply(sinb[numpy.newaxis,numpy.newaxis,:], sinalpha*sigmaphi) )**2 + ( numpy.multiply(cosb[numpy.newaxis,numpy.newaxis,:], sigmaz) )**2
    sigmaappl2 = numpy.divide(sigmal2, d2[:,numpy.newaxis,numpy.newaxis])
    sigmaappb2 = numpy.divide(sigmab2, d2[:,numpy.newaxis,numpy.newaxis])
#    sigmaappl2 = ? # vestigial
#    sigmaappb2 = ? # vestigial
    return (vrelappl, sigmaappl2, sigmaappb2)
